function [data] = osl_reduce_data_to_visualize(S)
%
% Takes a full erf [voxels x time] dataset or a tf [voxels x time x freq]
% data set and depending on the flags in 'params':
% - reduces the number of samples in the time domain by downsampling (so it
% can be quickly exported to fsl-view, for example)
% - can either select a frequency band, a range of frequency bands, or
% average over a few frequency bands
%
% Some input arguments (all optional; see script for defaults)
%
% FREQ
% - S.average_freqs (1/0): average together different frequency bands?
% - S.freq_bins (e.g. [3 4 5]). If averaging freqs, these will be the ones
%   averaged over.  If S.average_freqs = 0, just selects these freqs. The
%   default if S.freq_bins is unset is to select all bins.
%
% TIME
% - S.reduce_time (1/0): resample the data in time?
%   -> S.resample_method:
%       -> 'fft_ds'   : downsample using fft method in osl_spm_resample
%       -> 'collapse' : average together all times over the specified time
%       range
%   -> S.time_ds : downsampling factor in time.
% - S.time_range (e.g. [0.35 1.2]).  Cut down the data to this time range
%   only.  NB this dtermines the 'collapsing window' for
%   "S.resample_method  = 'collapse'".  Default if unset is to use all
%   timepoints.
%
% SPACE
% - S.doVoxelAvg (1/0): average over voxels?
% - S.voxels (e.g. [213:467]) : a list of the indices of the voxels over
%   which to average. If unset, default is to average over all of the
%   voxels in the input array.

data=S.data;
params=S;

% FLAGS -------------------------------------------------------------------
try tmp = params.average_freqs; catch, params.average_freqs = 0; end
try tmp = params.reduce_time; catch, params.reduce_time = 0; end
try tmp = params.resample_method; catch, params.resample_method = 'fft_ds'; end

% which data to use?
try tmp = params.freq_bins; catch, params.freq_bins = 1:size(data.cope,4); end
try tmp = params.time_ds; catch, params.time_ds = 1; end
try tmp = params.time_range; allTime = 0; catch, allTime = 1; params.time_range = [S.data.times(1) S.data.times(end)]; end

% average over certain voxels?  Default is to include all voxels (could be
% passing in data that was already masked, so this may sometimes make
% sense)
try tmp = params.doVoxelAvg; catch, params.doVoxelAvg = 0; end
try tmp = params.voxels; voxelsUnspecified = 0; catch, voxelsUnspecified = 1; params.voxels = 1:size(data.cope,1); end

try tmp = data.stdcope; is_stdcope=1; catch, is_stdcope=0; end

try, resolve_erf_sign_ambiguity=S.resolve_erf_sign_ambiguity; catch, Sin.resolve_erf_sign_ambiguity=0; end;

if size(data.cope,4) >1
    isTF = true;
else
    isTF = false;
end
% ------------------------------------------------------------------- FLAGS


% VOXELS ------------------------------------------------------------------
% if averaging over only certain voxels (to produce tf plots) do that now
if params.doVoxelAvg
    % fixed effect average over the specified voxels
    
    if(is_stdcope)
        %cope_var = data.stdcope(params.voxels,:,:,:).^2;
        %mean_cope_var = mean(cope_var,1);
        %data.stdcope = sqrt(mean_cope_var);
         
        dims=size(data.cope);
        dims(1)=1;
        datanew.cope=zeros(dims);
        datanew.stdcope=zeros(dims);
        
        
        if(Sin.resolve_erf_sign_ambiguity)
            if(size(data.cope,4)>1),
                error('resolve_erf_sign_ambiguity not compatible with TF data');
            end;

            for c=1:size(data.stdcope,3),
                for gc=1:size(data.stdcope,5),
                    tcs=data.cope(vox_coord,:,c,1,gc);
                    [a tc_pca]= pca(tcs,1);

                    if mean(tc_pca)<0,
                        tc_pca=-tc_pca;
                    end;

                    for ii=1:length(vox_coord),
                        ccs=corrcoef(tcs(ii,:),tc_pca);
                        if(ccs(1,2)<0),
                            tcs(ii,:)=-tcs(ii,:);
                        end;
                    end;    
                    data.cope(vox_coord,:,c,1,gc)=tcs;
                end;
            end;
        end;
        
        do_fixed_effects=1;
        for t=1:size(data.cope,2),
        for c=1:size(data.cope,3),
        for f=1:size(data.cope,4),
        for gc=1:size(data.cope,5),
            S=data.stdcope(params.voxels,t,c,f,gc).^2;
            z=ones(length(S),1);
            y=data.cope(params.voxels,t,c,f,gc);
            [gam, ~, covgam] = flame1(y,z,S,do_fixed_effects);
            datanew.cope(1,t,c,f,gc)=gam;
            datanew.stdcope(1,t,c,f,gc)=sqrt(covgam);
        end;
        end;
        end;
        end;

        data.cope=datanew.cope;
        data.stdcope=datanew.stdcope;
    else
        data.cope = mean(data.cope(params.voxels,:,:,:),1);    
    end;
    
    if voxelsUnspecified
        disp(['Averaging over all voxels (no voxel list given)...'])
    else
        disp(['Averaging over specified voxels...'])
    end
end % if params.doVoxelAvg
% ------------------------------------------------------------------ VOXELS




% TIME --------------------------------------------------------------------
% if downsampling in time, do that now
if params.reduce_time
    switch params.resample_method
        case 'fft_ds' % resample using the borrowed spm resample function (included as a subfunction)
            % NB may introduce wierd biases and edge effects, especially
            % with aggressive downsampling.  Use original data for stats!
            
            data.cope = osl_spm_resample(data.cope,params.time_ds);
            if(is_stdcope)
                cope_var = data.stdcope.^2;
                cope_var = osl_spm_resample(cope_var,params.time_ds);
                
                data.stdcope = sqrt(cope_var);
            end;
            data.times = linspace(min(data.times),max(data.times),size(data.cope,2));
            temp_res = 1/(data.times(2) - data.times(1));
            disp(['Downsampling in time: data now has ' num2str(size(data.cope,2)) ' time points.']);
            disp(['Temporal resolution: ' num2str(temp_res) ' Hz.'])
            time_indices=data.times>=params.time_range(1) & data.times<=params.time_range(2);
            data.cope=data.cope(:,time_indices,:,:,:);
            if(is_stdcope)
                data.stdcope=data.stdcope(:,time_indices,:,:,:);
            end;
            data.times=data.times(time_indices);
        case 'collapse' % collapse across the specified temporal interval
            time_indices=data.times>=params.time_range(1) & data.times<=params.time_range(2);
            data.times=data.times(time_indices);
            data.times2 = data.times;
            data.cope = data.cope(:,time_indices,:,:,:);
            data.cope    = mean(data.cope,2);
            if(is_stdcope)
                data.stdcope=data.stdcope(:,time_indices,:,:,:);
                cope_var     = data.stdcope.^2;
                cope_var     = mean(cope_var,2);
                data.stdcope = sqrt(cope_var);
            end;
            data.times   = mean(data.times);
            if allTime
                disp('Averaging across all time points');
            else
                disp(['Averaging across the time range ' num2str(data.times2(1)) ' to ' num2str(data.times2(end)) '.'])
            end
        otherwise
            error('''params.resample_method'' isn''t valid. Have a dig around in the function to see valid options.')
    end
else
    % just select the time bins requested
    time_indices=data.times>=params.time_range(1) & data.times<=params.time_range(2);
    data.times=data.times(time_indices);
    data.cope = data.cope(:,time_indices,:,:,:);
    if(is_stdcope)
        data.stdcope=data.stdcope(:,time_indices,:,:,:);
    end;
    disp(['Selected times ' num2str(data.times(1)) ' to ' num2str(data.times(end)) '.']);
end % if params_reduce_time
% -------------------------------------------------------------------- TIME



% FREQUENCY ---------------------------------------------------------------
% if averaging over frequency bins, do the average now
if isTF
    if params.average_freqs % to average over a few frequency bands
        data.cope = mean(data.cope(:,:,:,params.freq_bins),4);
        if(is_stdcope)
            cope_var = data.stdcope(:,:,:,params.freq_bins).^2;
            mean_cope_var = mean(cope_var,4);
            data.stdcope = sqrt(mean_cope_var);
        end;
        disp(['Averaging over ' num2str(data.frequencies(params.freq_bins)) ' Hz.'])
        data.frequencies = mean(data.frequencies(params.freq_bins));
    else % if we just want to select certain frequencies
        data.cope = data.cope(:,:,:,params.freq_bins);
        if(is_stdcope)
            data.stdcope = data.stdcope(:,:,:,params.freq_bins);
        end;
        disp(['Using ' num2str(data.frequencies(params.freq_bins)) ' Hz.'])
        data.frequencies = data.frequencies(params.freq_bins);
    end
else
    disp('Data is an erf, or single-frequency tf.')
    % if our data is an erf, or a tf but already compressed in the
    % frequency domain, then no need to subselect or average across
    % frequencies
end % if isTF
% --------------------------------------------------------------- FREQUENCY



% -------------------------------------------------------------------------
function resampled_signal = osl_spm_resample(signal,ds_factor)
% this is to downsample the time domain signal.  It's the spm
% function 'spm_resample'

N0          = size(signal,2);
N           = floor(N0/ds_factor);
ds_factor   = N/N0;
Y           = fftshift(fft(signal,[],2),2);
sy          = size(Y,2);
middle      = floor(sy./2)+1;

if ds_factor>1 % upsample
    N2 = floor((N-N0)./2);
    if N0/2 == floor(N0/2)
        Y(:,1) = []; % throw away non symmetric DFT coef
    end
    Y  = [zeros(size(Y,1),N2),Y,zeros(size(Y,1),N2)];
else % downsample
    N2 = floor(N./2);
    Y  = Y(:,middle-N2:middle+N2,:,:);
end
resampled_signal = ds_factor*ifft(ifftshift(Y,2),[],2);
